% Thus function solves the DSGE model under perfect foresight by solving a
% fixed point problem using the Newton-Raphson routine where the length of 
% the step size in each iteration is computed optimally using a simple grid
% search.
function [Z,f,diff,iter,errorMes] = FixedPointSolver_deltaOpt(Z0,x_t,y_tN,N,setupEPer)

% Unfold setupEPer
ny       = setupEPer.ny;
mx       = setupEPer.mx;
myx      = setupEPer.myx;
dispOn   = setupEPer.dispOn;

gridPoints = [0.0001 0.0001 0.001 0.01 0.1 0.2 0.4 0.6 0.8 1.0];
numGrid    = size(gridPoints,2);
% nx1+ny
nzdim2     = ny+mx;
Zvec       = reshape(Z0,N*nzdim2,1);
ZvecTmp    = zeros(N*nzdim2,numGrid);
fTmp       = zeros(N*nzdim2,numGrid);
diffTmp    = zeros(numGrid,1);
f          = DSGEforesight_N(Z0,x_t,y_tN,N,setupEPer);
Q          = sqrt(sum(real(f).^2)+sum(imag(f).^2));
Qold       = Q*100;
invW       = zeros(ny+mx,ny+mx,N);
fx11       = zeros(ny+mx,mx,N);
fx12       = zeros(ny+mx,myx,N);
fyp        = zeros(ny+mx,ny,N);
iter       = 0;
errorMes   = 0;
diff       = 1;
while abs(diff) > setupEPer.tolf && iter < setupEPer.MaxIter && ~isnan(diff) && Q < Qold
    Qold = Q;
    iter = iter + 1;
    % The gradient
    if setupEPer.JacobianOption == 1
        epsValue = 1e-8;
        dimf     = N*(ny+mx);
        fp       = zeros(dimf,dimf);
        fm       = zeros(dimf,dimf);
        for i=1:dimf
            % Positive step
            ZvecPlus      = Zvec;
            ZvecPlus(i,1) = Zvec(i,1) + epsValue;
            fp(:,i)       = DSGEforesight_N(reshape(ZvecPlus,ny+mx,N),x_t,y_tN,N,setupEPer);
            % Negative step
            ZvecMinus      = Zvec;
            ZvecMinus(i,1) = Zvec(i,1) - epsValue;
            fm(:,i)        = DSGEforesight_N(reshape(ZvecMinus,ny+mx,N),x_t,y_tN,N,setupEPer);
        end
        % Gradient
        J = (fp - fm)/(2*epsValue);
        delta  = J\(-f);
        % For debugging
        %Janal = analJacobian(reshape(Zvec,ny+mx,N),x_t,y_tN,N,setupEPer);;
        %max(max(abs(Janal-J)))
    elseif setupEPer.JacobianOption == 2
        J = analJacobian(reshape(Zvec,ny+mx,N),x_t,y_tN,N,setupEPer);
        if condest(J) < 1e12
            delta  = J\(-f);
        else
            delta  = zeros(N*(ny+mx),1);
        end
    elseif setupEPer.JacobianOption == 3
        updateJ = 1;
        [delta,invW,fx11,fx12,fyp] = nextStepFixedPointSolverBoucekkine(reshape(Zvec,ny+mx,N),x_t,y_tN,N,...
            reshape(f,ny+mx,N),invW,fx11,fx12,fyp,updateJ,setupEPer);
    end
    
    %The grid search
    for i=1:numGrid
        ZvecTmp(:,i) = Zvec + real(delta)*gridPoints(1,i);
        fTmp(:,i)    = DSGEforesight_N(reshape(ZvecTmp(:,i),nzdim2,N),x_t,y_tN,N,setupEPer);
        diffTmp(i,1) = sqrt(sum(real(fTmp(:,i)).^2+imag(fTmp(:,i)).^2));
        if isnan(diffTmp(i,1)) 
            diffTmp(i,1) = 1e60;
        end
    end
    [~,idx] = sort(diffTmp);
    f       = fTmp(:,idx(1,1));
    Q       = diffTmp(idx(1,1),1);
    scaling = gridPoints(1,idx(1,1));
    ZvecNew = ZvecTmp(:,idx(1,1));
    diff    = max(abs(ZvecNew - Zvec));
    
    % Updating
    if ~isnan(diff) && Q < Qold
        Zvec    = ZvecNew;
    end
     % We stop if the residuals of f0 are smaller than residualMax
    if setupEPer.residualMax > 0
        if max(abs(f(1:ny+mx,1))) < setupEPer.residualMax
            diff = 0;
        end
    end
    if dispOn == 1
        disp(['Q                = ', num2str(Q)])
        disp(['dQ               = ', num2str(Q-Qold)])
        disp(['Change in Z      = ', num2str(diff)])
        disp(['Scaling of delta = ', num2str(scaling)]);
        disp(['Iteration        = ', num2str(iter)]);
        disp(' ');
    end
end
if setupEPer.residualMax > 0
    % For the hybrid case where the criteria implied by setupEPer.residualMax is sufficient.
    if abs(diff) > setupEPer.tolf || iter > setupEPer.MaxIter || any(isnan(diff))
        errorMes = 1;
    end
else
    % Not for the hybrid case
    if abs(diff) > setupEPer.tolf || iter > setupEPer.MaxIter || any(isnan(diff)) || max(abs(f)) > setupEPer.tolf
        errorMes = 1;
    end
end
Z = reshape(Zvec,nzdim2,N);
end

